### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
##  This program conducts the Fama-French return regressions for     ##
##  The Business Cycle Dynamics of the Wealth Distribution           ##
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###


# Preliminaries -----------------------------------------------------------

rm(list=ls())

# Set directories ---------------------------------------------------------

#m1 <- Sys.getenv("LOGNAME")

home_directory = "[your home directory]"
project_directory <- file.path(home_directory, "replication_table2")
master_scripts_directory<-file.path(home_directory, "code") 

# Load in data ------------------------------------------------------------

library(stargazer)

dfa_returns <- read_csv(file.path(project_directory, "data/dfa_returns_for_FF.csv")) %>% 
  pivot_wider(names_from = cat, values_from = r_p) %>% 
  mutate(date = as.Date(as.yearqtr(date), frac=1)) %>% 
  mutate_if(is.numeric, ~.*100)

pvt_eqty <- read_csv(file.path(project_directory, "data/splpeqty_q.csv")) 
peq_date = pvt_eqty['date']
peq = pvt_eqty['PX_LAST']
colnames(peq) <- c("peq")
peq_all = data.frame(peq_date,peq)

smb_hml <- read_csv(file.path(project_directory, "data/ff_smb_hml_m.csv")) 
smb_hml_date = smb_hml['date']
smb = smb_hml['SMB']
hml = smb_hml['HML']
colnames(smb) <- c("smb")
colnames(hml) <- c("hml")
smb_hml_all = data.frame(smb_hml_date,smb,hml)

library(httr)    
set_config(use_proxy(url="198.35.160.54",port=8080))

library(fredr)
# replace this with your API key 
fredr_set_key("your API key")

fred_series <- list(
  series_id = c("WGS1YR", # one-year treasury yield (risk free rate)
                "BAMLCC0A0CMTRIV", # corporate bond index
                "DGS10" # ten year treasury rate 
  )
) 

fred <- purrr::pmap_dfr(
  .l = fred_series, 
  .f = ~fredr(series_id = .x, frequency = "q", observation_start = as.Date("1989-06-30"))
) %>% 
  select(date, series_id, value) %>% 
  pivot_wider(names_from = series_id) %>% 
  arrange(date) %>% 
  dplyr::rename(
    riskfree = WGS1YR,
    corpbond = BAMLCC0A0CMTRIV,
    tenyear = DGS10
  ) %>% 
  mutate(date = as.Date(as.yearqtr(date),frac=1))

fof <- fame2df(
  c(corelogic = "FI075035243.Q"),
  db = "fof",
  start = 19890630
)

us_q <-  fame2df(
  c(
    savrate = "SVPC%YPD.Q",
    unemp = "RUC.Q"
  ),
  db = "us",
  start = 19890630
)

fred_monthly_series <- list(
  series_id = c("FEDFUNDS",
                "SP500",
                "CEU0500000003"
  )
) 

us_m <- fame2df(
  serlist = c(
    fedfunds = "RIFSPFF_N.M",
    snp500 = "PJSPCC_N.M",
    avgearnings = "WRTH_I12Y@SGV_N.M" 
  ),
  db = "us",
  start = 19840930
) %>% 
  mutate(quarter = as.yearqtr(date)) %>% 
  group_by(quarter) %>%
  summarise_if(is.numeric, last) %>% 
  mutate(date = as.Date(quarter, frac=1)) %>% 
  select(date, where(is.numeric))

quart_return <- function(x, n=1) {
  (x - dplyr::lag(x, n)) / dplyr::lag(x, n) * 100
}

regdata <- Reduce(full_join, list(fred, fof, us_q, us_m)) 
regdata <- regdata %>% ungroup() %>% filter(date <= "2023-09-30")
regdata <- regdata %>% 
  arrange(date) %>%
  transmute(
    date = date,
    riskfree = ((1 + riskfree/100)^(1/4)-1) * 100,
    tenyear = ((1 + tenyear/100)^(1/4)-1) * 100,
    tenyear_5yrchange = tenyear - dplyr::lag(tenyear, 20),
    snp500 = quart_return(snp500) - riskfree,
    pvteqty = quart_return(peq) - riskfree,
    #smb_ff = quart_return(smb) - riskfree,
    smb_ff = smb - riskfree,
    hml_ff = hml - riskfree,
    corelogic = quart_return(corelogic) - riskfree,
    unemp = unemp - dplyr::lag(unemp),
    fedfunds_5yrchange = fedfunds - dplyr::lag(fedfunds, 20),
    fedfunds_1qtrchange = fedfunds - dplyr::lag(fedfunds),
    savrate = savrate - dplyr::lag(savrate),
    wagegrowth = quart_return(avgearnings),
    corpbond = quart_return(corpbond) - riskfree
  ) %>% 
  filter(date >= as.Date("1989-12-31")) %>% 
  full_join(dfa_returns, by = "date") %>% 
  mutate(
    TopPt1 = TopPt1 - riskfree,
    pct99to999 = pct99to999 - riskfree,
    pct90to99 = pct90to99 - riskfree,
    pct70to90 = pct70to90 - riskfree,
    pct50to70 = pct50to70 - riskfree,
    Bottom50 = Bottom50 - riskfree
  ) 

dfa_groups <- names(dfa_returns %>% select(where(is.numeric)))

regset <- list()

# Top 0.1% 

model1_toppt1 <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange, data=regdata)
model2_toppt1 <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange, data=regdata)
model3_toppt1 <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond, data=regdata)

model4_toppt1 <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + savrate, data=regdata)
model5_toppt1 <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + unemp, data=regdata)
model6_toppt1 <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + wagegrowth, data=regdata)
model7_toppt1 <- lm(TopPt1 ~ snp500 + corelogic + corpbond, data=regdata)
model8_toppt1 <- lm(TopPt1 ~ snp500 + corelogic + tenyear, data=regdata)
model9_toppt1 <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + tenyear, data=regdata)


model10_toppt1 <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond + savrate + unemp + wagegrowth + tenyear, data=regdata)

# 99 - 99.9% 

model1_pct99to999 <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange, data=regdata)
model2_pct99to999 <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange, data=regdata)
model3_pct99to999 <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond, data=regdata)
model4_pct99to999 <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + savrate, data=regdata)
model5_pct99to999 <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + unemp, data=regdata)
model6_pct99to999 <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + wagegrowth, data=regdata)
model7_pct99to999 <- lm(pct99to999 ~ snp500 + corelogic + corpbond, data=regdata)
model8_pct99to999 <- lm(pct99to999 ~ snp500 + corelogic + tenyear, data=regdata)
model9_pct99to999 <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + tenyear, data=regdata)

model10_pct99to999<- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond + savrate + unemp + wagegrowth + tenyear, data=regdata)

# 90 - 99% 

model1_pct90to99 <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange, data=regdata)
model2_pct90to99 <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange, data=regdata)
model3_pct90to99 <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond, data=regdata)
model4_pct90to99 <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + savrate, data=regdata)
model5_pct90to99 <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + unemp, data=regdata)
model6_pct90to99 <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + wagegrowth, data=regdata)
model7_pct90to99 <- lm(pct90to99 ~ snp500 + corelogic + corpbond, data=regdata)
model8_pct90to99 <- lm(pct90to99 ~ snp500 + corelogic + tenyear, data=regdata)
model9_pct90to99 <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + tenyear, data=regdata)

model10_pct90to99 <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond + savrate + unemp + wagegrowth + tenyear, data=regdata)

# 70 - 90% 

model1_pct70to90 <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange, data=regdata)
model2_pct70to90 <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange, data=regdata)
model3_pct70to90 <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond, data=regdata)
model4_pct70to90 <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + savrate, data=regdata)
model5_pct70to90 <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + unemp, data=regdata)
model6_pct70to90 <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + wagegrowth, data=regdata)
model7_pct70to90 <- lm(pct70to90 ~ snp500 + corelogic + corpbond, data=regdata)
model8_pct70to90 <- lm(pct70to90 ~ snp500 + corelogic + tenyear, data=regdata)
model9_pct70to90 <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + tenyear, data=regdata)

model10_pct70to90 <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond + savrate + unemp + wagegrowth + tenyear, data=regdata)

# 50 - 70% 

model1_pct50to70 <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange, data=regdata)
model2_pct50to70 <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange, data=regdata)
model3_pct50to70 <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond, data=regdata)
model4_pct50to70 <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + savrate, data=regdata)
model5_pct50to70 <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + unemp, data=regdata)
model6_pct50to70 <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + wagegrowth, data=regdata)
model7_pct50to70 <- lm(pct50to70 ~ snp500 + corelogic + corpbond, data=regdata)
model8_pct50to70 <- lm(pct50to70 ~ snp500 + corelogic + tenyear, data=regdata)
model9_pct50to70 <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + tenyear, data=regdata)

model10_pct50to70 <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond + savrate + unemp + wagegrowth + tenyear, data=regdata)


# Bottom 50%  

model1_bottom50 <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange, data=regdata)
model2_bottom50 <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange, data=regdata)
model3_bottom50 <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond, data=regdata)
model4_bottom50 <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + savrate, data=regdata)
model5_bottom50 <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + unemp, data=regdata)
model6_bottom50 <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + wagegrowth, data=regdata)
model7_bottom50 <- lm(Bottom50 ~ snp500 + corelogic + corpbond, data=regdata)
model8_bottom50 <- lm(Bottom50 ~ snp500 + corelogic + tenyear, data=regdata)
model9_bottom50 <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + tenyear, data=regdata)


model10_bottom50 <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + fedfunds_5yrchange + corpbond + savrate + unemp + wagegrowth + tenyear, data=regdata)

# Comment out -- we end up using model 5 of the pre-COVID period (found below)
#sink(file.path(project_directory, "regtables_mar23.txt"))
#stargazer(model1_toppt1, model2_toppt1, model3_toppt1, model4_toppt1, model5_toppt1, model6_toppt1, model7_toppt1, model8_toppt1, model9_toppt1, model10_toppt1, df = FALSE)
#stargazer(model1_pct99to999, model2_pct99to999, model3_pct99to999, model4_pct99to999, model5_pct99to999, model6_pct99to999, model7_pct99to999, model8_pct99to999, model9_pct99to999, model10_pct99to999, df = FALSE)
#stargazer(model1_pct90to99, model2_pct90to99, model3_pct90to99, model4_pct90to99, model5_pct90to99, model6_pct90to99, model7_pct90to99, model8_pct90to99, model9_pct90to99, model10_pct90to99, df = FALSE)
#stargazer(model1_pct70to90, model2_pct70to90, model3_pct70to90, model4_pct70to90, model5_pct70to90, model6_pct70to90,model7_pct70to90, model8_pct70to90, model9_pct70to90, model10_pct70to90 ,df = FALSE)
#stargazer(model1_pct50to70, model2_pct50to70, model3_pct50to70, model4_pct50to70, model5_pct50to70, model6_pct50to70, model7_pct50to70, model8_pct50to70, model9_pct50to70, model10_pct50to70 ,df = FALSE)
#stargazer(model1_bottom50, model2_bottom50, model3_bottom50, model4_bottom50, model5_bottom50, model6_bottom50, model7_bottom50, model8_bottom50, model9_bottom50, model10_bottom50, df = FALSE)
#sink()

#write_csv(regdata, file.path(project_directory, "regdata_mar23.csv"))

## final Latex table for paper 

regdata_precovid <- regdata %>% 
  filter(date <= as.Date("2019-12-31"))

reg_toppt1 <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata)
reg_pct99to999 <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata)
reg_pct90to99 <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata)
reg_pct70to90<- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata)
reg_pct50to70 <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata)
reg_bottom50 <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata)

reg_toppt1_precovid <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_precovid)
reg_pct99to999_precovid <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_precovid)
reg_pct90to99_precovid <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_precovid)
reg_pct70to90_precovid <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_precovid)
reg_pct50to70_precovid <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_precovid)
reg_bottom50_precovid <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_precovid)


sink(file.path(project_directory, "regtables_model5_mar23.txt"))
stargazer(reg_toppt1_precovid, reg_pct99to999_precovid, reg_pct90to99_precovid, 
          title = "Pre-COVID-19 Fama-French Regression", font.size = "footnotesize")
stargazer(reg_pct70to90_precovid, reg_pct50to70_precovid, reg_bottom50_precovid,
          title = "Pre-COVID-19 Fama-French Regression", font.size = "footnotesize")
stargazer(reg_toppt1, reg_pct99to999, reg_pct90to99, 
          title = "Full sample (including COVID-19) Fama French regression", font.size = "footnotesize")
stargazer(reg_pct70to90, reg_pct50to70, reg_bottom50,
          title = "Full sample (including COVID-19) Fama French regression", font.size = "footnotesize")
sink()

###### Why is R23 high in 99-99.9 group? S&P
reg_pct99to999_precovid <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_precovid)
reg_pct99to999_pc_sp <- lm(pct99to999 ~ snp500 , data=regdata_precovid)
reg_pct99to999_pc_addCL <- lm(pct99to999 ~ snp500 + corelogic , data=regdata_precovid)
reg_pct99to999_pc_addFF <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange , data=regdata_precovid)
reg_pct99to999_pc_add10 <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange , data=regdata_precovid)
reg_pct99to999_pc_CL <- lm(pct99to999 ~ corelogic , data=regdata_precovid)
reg_pct99to999_pc_FF <- lm(pct99to999 ~ fedfunds_1qtrchange , data=regdata_precovid)
reg_pct99to999_pc_10 <- lm(pct99to999 ~ tenyear_5yrchange , data=regdata_precovid)
reg_pct99to999_pc_ur <- lm(pct99to999 ~ unemp , data=regdata_precovid)

summary(reg_pct99to999_precovid)$r.squared
summary(reg_pct99to999_pc_sp)$r.squared
summary(reg_pct99to999_pc_addCL)$r.squared
summary(reg_pct99to999_pc_addFF)$r.squared
summary(reg_pct99to999_pc_add10)$r.squared

summary(reg_pct99to999_pc_CL)$r.squared
summary(reg_pct99to999_pc_FF)$r.squared
summary(reg_pct99to999_pc_10)$r.squared
summary(reg_pct99to999_pc_ur)$r.squared

#stargazer(reg_toppt1_precovid, reg_pct99to999_precovid, reg_pct90to99_precovid, reg_pct70to90_precovid, reg_pct50to70_precovid, reg_bottom50_precovid,
#          title = "Pre-COVID-19 Fama-French Regression", font.size = "footnotesize")
#stargazer(reg_toppt1, reg_pct99to999, reg_pct90to99, reg_pct70to90, reg_pct50to70, reg_bottom50,
#          title = "Full sample (including COVID-19) Fama French regression", font.size = "footnotesize")
#sink()

reg_toppt1_pcvd_pe <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_precovid)
reg_pct99to999_pcvd_pe <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_precovid)
reg_pct90to99_pcvd_pe <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_precovid)
reg_pct70to90_pcvd_pe <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_precovid)
reg_pct50to70_pcvd_pe <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_precovid)
reg_bottom50_pcvd_pe <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_precovid)

sink(file.path(project_directory, "regtables_model5_PE_mar23.txt"))
stargazer(reg_toppt1_pcvd_pe, reg_pct99to999_pcvd_pe, reg_pct90to99_pcvd_pe,
          title = "Pre-COVID-19 Fama-French Regression (w/ Pvt. Eqty)", font.size = "footnotesize")
stargazer(reg_pct70to90_pcvd_pe, reg_pct50to70_pcvd_pe, reg_bottom50_pcvd_pe,
          title = "Pre-COVID-19 Fama-French Regression (w/ Pvt. Eqty)", font.size = "footnotesize")
sink()

#stargazer(reg_toppt1_pcvd_pe, reg_pct99to999_pcvd_pe, reg_pct90to99_pcvd_pe, reg_pct70to90_pcvd_pe, reg_pct50to70_pcvd_pe, reg_bottom50_pcvd_pe,
#          title = "Pre-COVID-19 Fama-French Regression", font.size = "footnotesize")
#stargazer(reg_toppt1, reg_pct99to999, reg_pct90to99, reg_pct70to90, reg_pct50to70, reg_bottom50,
#          title = "Full sample (including COVID-19) Fama French regression", font.size = "footnotesize")
#sink()

# Pre-Covid, 2004 onwards (2003 is when pvt equity series begins)
regdata_pc_pe <- regdata_precovid %>% 
  filter(date >= as.Date("2004-03-31"))

# 2003-2019 period -- repeat regression w/o pvt equity
reg_toppt1_pc_nope <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_pc_pe)
reg_pct99to999_pc_nope <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_pc_pe)
reg_pct90to99_pc_nope <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_pc_pe)
reg_pct70to90_pc_nope <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_pc_pe)
reg_pct50to70_pc_nope <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_pc_pe)
reg_bottom50_pc_nope <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp, data=regdata_pc_pe)

# 2003-2019 period --  regression w/ pvt equity
reg_toppt1_pc_pe <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_pc_pe)
reg_pct99to999_pc_pe <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_pc_pe)
reg_pct90to99_pc_pe <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_pc_pe)
reg_pct70to90_pc_pe <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_pc_pe)
reg_pct50to70_pc_pe <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_pc_pe)
reg_bottom50_pc_pe <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + pvteqty$peq, data=regdata_pc_pe)

sink(file.path(project_directory, "regtables_model5_20042019pd.txt"))
stargazer(reg_toppt1_pc_pe, reg_pct99to999_pc_pe, reg_pct90to99_pc_pe,
          title = "2004-2019 Fama-French Regression (w/ Pvt. Eqty)", font.size = "footnotesize")
stargazer(reg_pct70to90_pc_pe, reg_pct50to70_pc_pe, reg_bottom50_pc_pe,
          title = "2004-2019 Fama-French Regression (w/ Pvt. Eqty)", font.size = "footnotesize")
stargazer(reg_toppt1_pc_nope, reg_pct99to999_pc_nope, reg_pct90to99_pc_nope,
          title = "2004-2019 Fama-French Regression (w/o Pvt. Eqty)", font.size = "footnotesize")
stargazer(reg_pct70to90_pc_nope, reg_pct50to70_pc_nope, reg_bottom50_pc_nope,
          title = "2004-2019 Fama-French Regression (w/o Pvt. Eqty)", font.size = "footnotesize")
sink()

#SMB and HML  
reg_toppt1_pc_smb <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + smb_ff$smb, data=regdata_precovid)
reg_pct99to999_pc_smb <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + smb_ff$smb, data=regdata_precovid)
reg_pct90to99_pc_smb <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + smb_ff$smb, data=regdata_precovid)
reg_pct70to90_pc_smb <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + smb_ff$smb, data=regdata_precovid)
reg_pct50to70_pc_smb <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + smb_ff$smb, data=regdata_precovid)
reg_bottom50_pc_smb <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + smb_ff$smb, data=regdata_precovid)

reg_toppt1_hml <- lm(TopPt1 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + hml_ff$hml, data=regdata_precovid)
reg_pct99to999_hml <- lm(pct99to999 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + hml_ff$hml, data=regdata_precovid)
reg_pct90to99_hml <- lm(pct90to99 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + hml_ff$hml, data=regdata_precovid)
reg_pct70to90_hml <- lm(pct70to90 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + hml_ff$hml, data=regdata_precovid)
reg_pct50to70_hml <- lm(pct50to70 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + hml_ff$hml, data=regdata_precovid)
reg_bottom50_hml <- lm(Bottom50 ~ snp500 + corelogic + fedfunds_1qtrchange + tenyear_5yrchange + unemp + hml_ff$hml, data=regdata_precovid)

